Subject-specific one-dimensional fluid dynamics model of chronic thromboembolic pulmonary hypertension

Chronic thromboembolic pulmonary hypertension (CTEPH) develops due to the accumulation of blood clots in the lung vasculature that obstruct flow and increase pressure. The mechanobiological factors that drive progression of CTEPH are not understood, in part because mechanical and hemodynamic changes in the pulmonary vasculature due to CTEPH are not easily measurable. Using previously published hemodynamic measurements and imaging from a large animal model of CTEPH, we developed a subject-specific one-dimensional (1D) computational fluid dynamic (CFD) models to investigate the impact of CTEPH on pulmonary artery stiffening, time averaged wall shear stress (TAWSS), and oscillatory shear index (OSI). Our results demonstrate that CTEPH increases pulmonary artery wall stiffness and decreases TAWSS in extralobar (main, right and left pulmonary arteries) and intralobar vessels. Moreover, CTEPH increases the percentage of the intralobar arterial network with both low TAWSS and high OSI. This subject-specific experimental-computational framework shows potential as a predictor of the impact of CTEPH on pulmonary arterial hemodynamics and pulmonary vascular mechanics. By leveraging advanced modeling techniques and calibrated model parameters, we predict spatial distributions of flow and pressure, from which we can compute potential physiomarkers of disease progression, including the combination of low mean wall shear stress with high oscillation. Ultimately, this approach can lead to more spatially targeted interventions that address the needs of individual CTEPH patients.


Introduction
Chronic thromboembolic pulmonary hypertension (CTEPH) occurs when blood clots lodge in the small blood vessels of the lung and fail to fully dissolve or resolve over time. The resulting blood ow obstruction leads to an elevation in blood pressure within the pulmonary circulation. If left untreated, CTEPH causes right ventricular failure and death. With early diagnosis, however, CTEPH can be successfully treated with surgery (Feinstein et al. 2001). CTEPH is diagnosed by invasive right heart catheterization (RHC) to measure mean pulmonary artery pressure (mPAP) ≥ 20 mmHg and pulmonary capillary wedge pressure (PCWP) ≤ 15 mmHg, with subsequent imaging to con rm clot burden (Lang et al. 2014).
A key knowledge gap in CTEPH is the drivers of clot growth and dissolution over time (Tsubata et al. 2023). Vascular mechanobiological factors are thought to play an important role. Abnormalities in the magnitude and directions of wall shear stress (WSS) contribute to pathological processes such as endothelial dysfunction, in ammation, and impaired vasodilation (Allen et al. 2023)WSS in small pulmonary arteries likely impact clot organization and growth as well as dissolution (Tsubata et al. 2023). Given that pulmonary artery stiffness increases with pulmonary hypertension development (Wang and Chesler 2011) and alters ow and thus shear stress dynamics (Su et al. 2013), pulmonary artery stiffness may also contribute to CTEPH progression. Therefore, computing WSS and vascular mechanics throughout the pulmonary arterial network can provide useful information regarding likely sites for disease progression.
Large animal models of CTEPH have advanced our understanding of the pathophysiology of CTEPH (Mulchrone et al. 2019;Stam et al. 2021), but typically do not report local mechanical stimuli for disease progression such as WSS. Image-based computational uid dynamics (CFD) modeling can compute pulmonary artery blood ow with high spatial resolution (Burrowes et al. 2017;Bordones et al. 2018; Wang et al. 2020), from which WSS can be calculated. The combination of experiments in large animal models of CTEPH with CFD analysis enables estimation of mechanobiological stimuli in sections of the pulmonary arterial network inaccessible to measurement. Moreover, this combination of techniques can be used to compute stimuli that cannot be measured from imaging alone, such as pressure distribution.
Here, we use an image-based nonlinear one-dimensional (1D) CFD model to predict hemodynamics throughout the pulmonary arterial network before and after CTEPH. The large animal data were previously published by (Mulchrone et al. 2019). We calibrate our CFD model to the hemodynamic and imaging data available and then use the results to quantify the impact of CTEPH on metrics of WSS including time-averaged WSS (TAWSS) and oscillatory shear index (OSI). We also quantify the impact of CTEPH on wall stiffness. This synergistic experimental-computational framework gives insight into the extralobar and intralobar mechanobiological effects of CTEPH with subject speci city and provides a future tool for understanding the local determinants of disease progression seen clinically.

Animal Study
All experimental data were obtained retrospectively; a detailed protocol and complete reporting of measurements can be found in (Mulchrone et al. 2019). Brie y, CTEPH was induced by injecting multiple microspheres into the pulmonary arteries (PAs) using an Page 3/17 indwelling catheter following methods established by Hori et al. (Hori et al. 2012) in ve adult male canines (12 ± 1 kg body weight).
Magnetic resonance imaging (MRI) and RHC were performed under anesthesia sustained by administering 1-3% iso urane in 100% oxygen, and ventilation regulated to maintain end-tidal CO 2

CFD Model Geometry
We created a 1D arterial network model as the computational domain for each subject pre-and post-CTEPH, shown for a representative subject in Fig. 1. First, from sagittal plane MR images (Fig. 1a), 3D segmentations of extralobar and intralobar PA geometries were created using 3D slicer (Fedorov et al. 2012) (Fig. 1b). We then used the vascular Modeling Toolkit to convert the 3D segmentations to 1D centerline networks (Fig. 1c) (Antiga et al. 2003). Finally, we post-processed the speci c networks using custom MATLAB (Natick, MA) software to create the 1D CFD mathematical domain for simulations (Colebank et al. 2021). To ensure agreement between the dynamic extralobar PA area in the model and the time-averaged extralobar PA area from MRI, we scaled the extralobar PA area in the network geometry at diastole to match the MRI data at diastole.

Governing Equations
We use time-dependent, 1D uid dynamics equations to simulate hemodynamics in the pulmonary arterial networks. We assume all blood vessels are straight, cylindrical, and impermeable. Blood ow through the arterial network is assumed to be axisymmetric, incompressible, Newtonian, and laminar such that the governing equations are: where (cm) and (s) represent the axial and temporal coordinates, and is the transmural blood pressure (mmHg) (Olufsen et al. 2000). is the volumetric ow (mL/s) and is cross sectional area (cm 2 ) where is the spatially and temporally-dependent radius of the artery (cm). The density ( ) and dynamic viscosity ( ) of blood are assumed to be constant and equal to 1.03 (g/mL) and 0.03 (Poise), respectively (Kamiya and Togawa 1980). Further we assume the uid velocity at the wall is equal to the velocity of the wall (no-slip condition).
The 1D theory requires an explicit expression for the uid velocity pro le. Here, we use the power-law pro le: where is the blood velocity (cm/s), is the mean velocity (cm/s), and is the power-law constant, which is set to 9 to provide a nearly uniform velocity pro le within the center of the vessel consistent with a high Womersley number (Alastruey et al. 2008; Spilker and Taylor 2010).

Wall Mechanics
Within the given 1D equations, there are three dependent variables, { . To close the system of equations, a constitutive model for the artery wall is required. We used a linear elastic model (Safaei et al. 2016), which assumes a linearly elastic, thin walled, isotropic, and incompressible material in the form of a cylinder. Under these circumstances, the constitutive law can be written as: 4 where is the circumferential Young's modulus (mmHg), (cm) is the wall thickness, and (cm 2 ) is the lumen area at the diastolic pressure (mmHg). The term accounts for both structural and load-dependent stiffening and is a potential physiomarker of disease severity. We calculated this term analytically for individual subjects pre-and post-CTEPH using systolic pressure and area in the MPA in Eq. (4). Finally, model equations were solved using two-step Lax-Wendroff method, detailed in prior work (Olufsen et al. 2000).

Boundary Conditions
Boundary conditions are speci ed at each vessel inlet and outlet. We used the as the inlet ow boundary condition in the MPA. At extralobar and intralobar vessel junctions, ow conservation and pressure continuity were enforced, where the subscript 'p' denotes the parent vessel that bifurcates into two daughter vessels 'd1' and 'd2'.
At each terminal intralobar vessel, we used 3 element Windkessel models to represent all vessels outside those obtained by segmentation. Each Windkessel model includes a proximal resistance , distal resistance , and compliance . Nominal where is the calibration pressure for , is computed pressure, is the calibration area, and is the computed area. For the time-series data terms, N is the length of the time vector spanning one cardiac cycle, denotes calibration ow data for , and is the computed time-series ow at time . We used the function lsqnonlin in MATLAB to identify the global minimum for our objective function based on ten randomized starting values (Qureshi et al. 2019). The calibration process is shown in Fig. 2.

Wall shear stress
The TAWSS and OSI were computed as where (6) and 7 throughout the arterial network. Here, (dyne/cm 2 is the spatially and temporally dependent arterial WSS and (s) is the cardiac cycle length. Note that whereas TAWSS and OSI in the large, extralobar arteries are straightforward to compare from baseline to CTEPH, the number and size of intrapulmonary arteries detectable in the imaging, and thus represented in the 1D model, varies from baseline and CTEPH. To enable comparison from baseline to CTEPH, we report TAWSS and OSI averaged over all intralobar arteries for each subject.

Statistical Test
Paired one-sided t-tests (in Excel) were used to assess signi cance with p < 0.05.

Results
The objective of this study was to formulate a comprehensive framework that leverages clinically measurable hemodynamic data and imaging to forecast non-clinically measurable hemodynamic and mechanical parameters relevant to pulmonary vascular disease progression. We used imaging data to build subject-speci c geometries of the pulmonary arterial network for uid dynamics simulation pre-and post-CTEPH induction. After calibrating to hemodynamics pre-and post-CTEPH, we performed 1D uid dynamics model simulations to investigate the impact of CTEPH on spatial and temporal distributions of wall shear stress and vascular mechanics.

Model calibration
Model ts to scalar calibration data { at baseline and CTEPH are shown in Fig. 3. Following the development of CTEPH, all subjects had increased systolic and diastolic pressures as expected ( Fig. 3a and b). Four out of ve subjects had increased systolic and diastolic areas. We attribute the one exception to low image resolution at baseline, which led to a non-signi cant group increase with CTEPH ( Fig. 3c and d).  Pulmonary arterial stiffness was analytically computed for each subject using { in Eq. (4). As shown in Fig. 5a, CTEPH signi cantly increased (p = 0.02), which can be attributed to the combination of structural and load-dependent stiffening.

Wall Shear Stress
With the model predictions of { , changes in TAWSS with CTEPH were computed (Fig. 6). In four out of ve subjects, TAWSS decreased at the MPA and RPA; in all subjects, TAWSS decreased in the LPA and most intralobar arteries. One subject exhibited an increase in MPA TAWSS (by 3%) and one subject exhibited an increase in RPA TAWSS (by 9%) following CTEPH. On average, there was a signi cant decrease in TAWSS in the MPA, LPA, RPA, and intralobar vessels (p = 0.02, p = 0.007, p = 0.04 and p = 0.02, respectively).
A three-dimensional map displaying the distribution of TAWSS across the entire arterial network for a representative subject at baseline and CTEPH is shown in Fig. 7. At baseline, TAWSS values ranged from 1.7 to 94.5 (dyne/cm 2 ) with a mean value of 22.6 whereas after CTEPH, TAWSS decreased, ranging from 1.4 to 64.3 (dyne/cm 2 ) with a mean value of 16.1.
Model ow simulations also enabled computation of OSI distribution. With CTEPH, OSI increased and decreased in the extralobar (MPA, LPA, and RPA) and intralobar arteries without a consistent pattern (Fig. 8). For several subjects, MPA and RPA OSI values were essentially zero.
To better understand the hemodynamic impact of CTEPH, we investigated changes in a combined metric using the spatial distribution of TAWSS and OSI. We identi ed regions of the network with both low TAWSS (< 5 dyne/cm 2 ) and high OSI (> 0.05); the TAWSS cutoff was based on previous literature (Li et al. 2009) and the OSI cutoff was based on the average OSI for all vessels across all subjects. As shown in Fig. 9, the percentage of regions that met this criterion, , signi cantly increased with CTEPH (p = 0.0004).

Discussion
In this study, we used an image-based nonlinear one-dimensional CFD model to predict arterial hemodynamic features that would be di cult, if not impossible, to measure experimentally and which may be clinically relevant to disease progression. Based on an existing set of imaging data pre-and post-CTEPH, we developed subject-speci c 1D arterial network models. We used an e cient multi target optimization scheme for calibration that incorporated experimentally measured MPA systolic and diastolic pressures and areas, as well as dynamic ow data at the LPA and RPA. Our approach yielded simulation results with excellent agreement with measured data, and predicted pulmonary arterial stiffness, TAWSS and OSI distribution pre-and post-CTEPH. The increase in pulmonary arterial stiffness post-CTEPH will drive upstream (i.e., right ventricular (Liu et al. 2022)) and downstream (i.e., pulmonary capillary (Wang and Chesler 2011)) remodeling to accelerate disease progression. Moreover, changes in WSS are relevant to the growth or dissolution of thrombi (Sakariassen et al. 2015). Thus, these model results allow for a more comprehensive understanding of the effects of CTEPH, enhancing our ability to study this condition and potentially improve patient care and treatment strategies.

Modeling and Optimization
Hemodynamics modeling is a powerful tool for investigating hypothesized mechanisms of pulmonary vascular disease development and progression. There has been limited application of these tools to CTEPH in comparison to other forms of pulmonary hypertension, terminal vessel in the out ow branches, were adjusted using in vivo measurements of pressure from RHC. Additionally, the linearized stiffness E of the large arteries was adjusted based on the relationship between pressure and diameter from MRI. Clipp  (time-series) calibration data, enhancing model accuracy and reliability. As shown in Fig. 4 for a representative animal and consistent with all animals, model predictions quantitatively matched measured pressures and ows (high 2 values).
As mentioned earlier, there are few studies that used computational modeling to predict hemodynamics in CTEPH. Colebank et al. (2021) developed 1D CFD models of pulmonary arterial networks based on computed tomography imaging in CTEPH patients. They simulated four different CTEPH disease cases to test physiological hypotheses related to remodeling. Tsubata et al. (2023) performed 3D CFD simulations to model subject-speci c data from healthy subjects, CTEPH patients and CTEPH patients treated with balloon pulmonary angioplasty. However, neither paper had an abundance of pre-and post-CTEPH data from individual subjects. A novel aspect of our study is that we used paired data sets pre-and post-CTEPH, which enables a subject-speci c interpretation of CTEPH progression.

Pulmonary Arterial Parameters
Pulmonary hypertension is known to increase pulmonary arterial stiffness, which is measured in clinical and preclinical studies using a wide range of direct and indirect methods. Su et al(2017) used wave intensity analysis to indirectly measure pulmonary arterial stiffness in 10 healthy subjects and 21 subjects with pulmonary hypertension (11 with PAH and 10 with CTEPH). They reported that the wave re ection index in these subjects was signi cantly higher than in healthy controls. We analytically determined the mechanical properties  At the MPA, we found that TAWSS signi cantly decreased following CTEPH in four of the ve subjects; in the one exception, that we attribute this to poor image quality at baseline that likely increased the estimated arterial network area. A signi cant reduction in LPA TAWSS was observed for all subjects in our study. The observed decrease in TAWSS can be attributed to a dramatic drop in ow to the left lung due to the severity of the obstruction. This nding is consistent with the suggestion that WSS reduction in the proximal PAs is an indicator of PAH disease severity (Yang et al. 2019). At the RPA, TAWSS decreased in all subjects except one (not the same subject for which MPA TAWSS increased). Overall, the decrease in TAWSS in the RPA was less dramatic and more variable than in the LPA (Fig. 6).
consistent with Pillalamarri et al(2021) who found that TAWSS increased with increasing distance from the heart based on a computational investigation of patient-speci c pulmonary hemodynamics across 32 subjects in distinct WHO groups of PH. The Pillalamarri et al. (Pillalamarri et al. 2021)study included one subject with CTEPH, in whom TAWSS decreased in the distal part of the network, consistent with our results.
The biological consequences of decreased WSS in the pulmonary circulation include altered endothelial cell signaling that drives vasoconstriction and the production of pro-in ammatory factors (Allen et al. 2023). Decreased TAWSS may also contribute to the synthesis and accumulation of collagen in the arterial wall and proliferation of smooth muscle cells (Ryu et al. 2022). This pulmonary arterial remodeling can result in a stiffer and less compliant pulmonary artery network (Friesen et al. 2019), which further exacerbates the severity of PH and can lead to right heart failure.

Oscillatory shear index
Like TAWSS, OSI can drive changes in endothelial cell structure and function that affect disease progression. In endothelial cells, high OSI is associated with increased reactive oxygen species, which can drive vasoconstriction and in ammation (Peiffer et al. 2013 In the intralobar arteries, there was no clear pattern regarding the impact of CTEPH on OSI (Fig. 8d). To further investigate temporal abnormalities in intralobar ow dynamics, we considered the combination of TAWSS and OSI (Fig. 9) represented by the dimensionless parameter . Tsubata et al. (2023) suggested that the combination of decreased WSS and increased OSI can contribute to thrombogenicity. Our results indicate that the percentage of vessels that had low TAWSS (< 5dyne/cm 2 ) and high OSI (> 0.05) increased from baseline to CTEPH in all subjects. When the vascular endothelium experiences a decrease in wall shear stress and an increase in OSI, it inhibits the production and release of nitric oxide, a potent vasodilator (Terada et al. 2016) As a result, the abnormal ow dynamics as seen in our study may be linked to impaired vasodilation and vasoconstriction in the pulmonary arterioles.

Limitations
Several limitations exist in our study, notably the small sample size, which contributed to heterogeneity. Limited image resolution in some subjects, especially one subject at baseline and a different subject after CTEPH, resulted in challenges extracting accurate subjectspeci c arterial networks and changes in those arterial networks with CTEPH. Since time-series pressure data were not recorded for either baseline or CTEPH conditions, we calibrated our model based on the scalar values of systolic and diastolic pressures (as well as scalar values of systolic and diastolic MPA area and time series LPA and RPA ow). The inclusion of time-series MPA pressure data would increase model delity. Finally, we assumed Blunt-like ow with a power law radial distribution for all arteries and held the power xed at g = 9. However, the range of Womersely number in the models was 2 to 40. As a consequence, in the small arteries, which will have lower Womersely numbers, TAWSS is overestimated. The impact on OSI is di cult to predict a priori; comparison to Tsubata et al. (2023) suggests OSI is underestimated in our model. Adapting g to local diameter or using an explicit Womersley pro le in in the 1D CFD model would increase model complexity and solution time but increase accuracy.

Conclusion
The impact of CTEPH on hemodynamic and vascular mechanical parameters that cannot be measured directly were examined using a subject-speci c 1D CFD modeling framework. We calibrated the model to RHC-and MRI-derived measurements of pressure, area, and ow from a previously published large animal model of CTEPH. Our analysis revealed that CTEPH increased arterial stiffness and decreased TAWSS, and that the combination of decreased TAWSS and increased OSI in intralobar arteries may correlate with CTEPH progression. Our experimental-computational framework shows potential as a patient-speci c simulator of both extralobar and intralobar pulmonary hemodynamics within the context of pulmonary hypertension. By leveraging advanced modeling techniques with model calibration, we can gain a more wholistic understanding of CTEPH, which could ultimately lead to more targeted interventions that address the speci c needs of individuals with CTEPH.  Schematic of model calibration process for Windkessel parameters with experimental pressure, area and ow data.  Measured hemodynamics and model predictions for a representative subject. Figure 3(a,b) shows model pressure predictions (solid lines) compared to measured systolic and diastolic pressure (dashed lines) for a representative animal in which mPAP increased from 20 to 51 mmHg from baseline (blue; left) to CTEPH (red; right).  ) show LPA ow and RPA ow, respectively, from measurements (black) and model predictions (blue, left at baseline and red, right with CTEPH). In this animal, a notable decrease in mean ow at LPA was observed after CTEPH induction, decreasing from 10.02 to 7.40 (mL/s), which represents a reduction of approximately 26%. Conversely, an increase in mean ow of approximately 46% was observed in the right lung. 3D map of TAWSS distribution across the entire arterial network for a representative animal at baseline (a) and CTEPH (b). This map demonstrates the changes in TAWSS that occurred in response to CTEPH induction throughout different segments of the network.